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ABSTRACT 

Motivated by a considerable scatter in the observationally inferred lifetimes of the embedded phase 
of star formation, we study the duration of the Class and Class I phases in upper-mass brown 
dwarfs and low-mass stars using numerical hydrodynamics simulations of the gravitational collapse of 
a large sample of cloud cores. We resolve the formation of a star/disk/envelope system and extend 
our numerical simulations to the late accretion phase when the envelope is nearly totally depleted of 
matter. We adopted a classification scheme of Andre et al. and calculate the lifetimes of the Class 
and Class I phases (rco and rci, respectively) based on the mass remaining in the envelope. When 
cloud cores with various rotation rates, masses, and sizes (but identical otherwise) are considered, our 
modeling reveals a sub-linear correlation between the Class lifetimes and stellar masses in the Class 
phase with the least-squares fit exponent m = 0.8 ± 0.05. The corresponding correlation between the 
Class I lifetimes and stellar masses in the Class I is super-linear with m = 1.2 ± 0.05. If a wider 
sample of cloud cores is considered, which includes possible variations in the initial gas temperature, 
cloud core truncation radii, density enhancement amplitudes, initial gas density and angular velocity 
profiles, and magnetic fields, then the corresponding exponents may decrease by as much as 0.3. The 
duration of the Class I phase is found to be longer than that of the Class phase in most models, with 
a mean ratio tqi/tco ~ 1.5-2. A notable exception are YSOs that form from cloud cores with large 
initial density enhancements, in which case rco may be greater than tqi. Moreover, the upper-mass 
(> 1.0 Mq) cloud cores with frozen-in magnetic fields and high cloud core rotation rates may have 
the tci/tco ratios as large as 3.0-4.0. We calculate the rate of mass accretion from the envelope onto 
the star/disk system and provide an approximation formula that can be used in semi- analytic models 
of cloud core collapse. 

Subject headings: circumstellar matter — planetary systems: protoplanetary disks — hydrodynamics 
— ISM: clouds — stars: formation 



1. INTRODUCTION 

Constraining the lifetimes associated with different 
phases of low-mass star formation has been tradition- 
ally one of the goals of research. In spite of much 
effort in this field, yet there is a considerable scat- 
ter among observationally estimated lifetimes of the 
embedded phase. Most previous observational stud- 
ies suggested the duration of the Class and Class I 
phase s to be within the 0.1 - 0.7 Myr range ( Wilking et al.1 



19891: iGreene et all Il994t iKenyon k Hartmannl 119951 : 
Visser et al.ll2002tlHatchell et al.ll2007l). with similar rel- 



ative lifetimes (j Visser et all 120021: lHatchell et all 125571 ), 
while some argue in favour of a significantly shorter 
lifetime for the Class phase, rco=0.01-0.06 Myr 
(jAndre k M ontmerlc 1994). In the most recent work 
involvi ng a large set of YSO s fro m different s t ar form ing 
clouds, lEnoch et all (|2009l ) and lEvans et all (|2009t ) ad- 
vocated for the lifetimes of 0.1-0.2 Myr for the Class 
sources and 0.44-0.54 Myr for the Class I ones. 

The aforementioned lifetimes may vary in part due 
to different estimation techniques and employed wave- 
lengths and in part due to observations being confined to 
different molecular clouds. The cornerstone assumption 
of steady star formation, adopted in most observational 
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studies, may break down if star formation relies heavily 
on external trigger ing, as may be the case in p Ophiuchi 
(| Visser et al.|[20 02j). Furthermore, the environment and 
local initial conditions determine how long the Class 
and Class I phases last. 

It is therefore interesting to compare the observation- 
ally inferred lifetimes with those predicted from numer- 
ical modeling of the c l oud co re collapse. For instance, 
iMasunaga k Inutsukal (|2000D performed radiation hy- 
drodynamics simulations of the gravitational collapse 
of spherically symmetric cloud cores and found, based 
on the synthetic spectral energy distribution, that the 
Class phase lasts around 2 x 10 4 yr. However, their 
modeling was limited to only two clo ud cores. Some 
work in this direction was also done by iFroebrich et al.1 
(2006), who employed smooth particle hydrodynamics 
simulations of the fragmentation and collapse of turbu- 
lent, self-gravitating clouds and found Class lifetimes 
to be of the order of (2 — 6) x 10 4 yr, significantly lower 
than those inferred from most observations. However, 
they did not provide estimates for the Class I lifetime, 
probably due to an enormous computational cost, and 
their lifetimes were inferred from the mass infall rates at 
~ 300 AU rather than from more customary diagnostics 
such as envelope masses or spectral properties. 

In this paper, we perform a comprehensive numer- 
ical study of the duration of the embedded phase of 
star formation in upper-mass brown dwarfs and low- 
mass stars. Using numerical hydrodynamics simulations, 
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we compute the gravitational collapse of a large sample 
of gravitationally unstable cloud cores with various ini- 
tial masses, rotation rates, gas temperatures, density en- 
hancement amplitudes, truncation radii, and strengths 
of frozen-in magnetic fields. The paper is organized as 
follows. In Sections [2] and [3l we briefly review the nu- 
merical model and initial conditions. In Section [U we 
discuss the adopted classification scheme of YSOs. In 
Section 15.11 we search for possible correlations between 
the Class and Class I lifetimes, from one hand, and 
stellar and cloud core masses, from the other hand. The 
effect of varying initial conditions is considered in Sec- 
tion [521 The model envelope depletion rates are tested 
against several empirical functions in Section [5] and the 
main results are summarized in Section [71 

2. MODEL DESCRIPTION 

Our numerical model is similar to that used recently to 
study the secular evolution of circumstellar disks and the 
mass accretion rates in T T auri stars and bro wn dwarfs 
(|Vorobvov fc Basul l2009aH bl; IVorobvovl l2009aft . For the 
reader's convenience, we briefly summarize the basic con- 
cept and equations. 

We make use of the thin-disk approximation to com- 
pute the long-term (~ 2 Myr) evolution of rotating, grav- 
itationally unstable cloud cores. We start our numer- 
ical integration in the pre-stellar phase, which is char- 
acterized by a collapsing starless cloud core, continue 
into the embedded phase, which sees the formation of 
a star/disk/envelope system, and terminate our simula- 
tions in the late accretion phase, when most of the cloud 
core has accreted onto the star/disk system. Once the 
disk has self-consistently formed, it occupies the inner- 
most regions of our numerical grid, while the inf ailing 
cloud core (the envelope) occupies the rest of the grid. 
This ensures that the mass accretion rate onto the disk 
M cnv is self-consistently determined by the dynamics of 
the gas in the envelope, rather than being simply intro- 
duced as a free parameter. 

We introduce a "sink cell" at r < 5 AU and impose 
a free inflow inner boundary condition. As the gravita- 
tional collapse of a cloud core ensues, the matter begins 
to freely flows through the sink cell. We monitor the gas 
surface density in the sink cell and when its value ex- 
ceeds a critical value S cr for the transition from isother- 
mal to adiabatic evolution (see Equation [5]), we intro- 
duce a central star. In the subsequent evolution, ninety 
per cent of the gas that crosses the inner boundary is 
assumed to land onto the central star plus the inner ax- 
isymmetric disk at r < 5 AU. The mass ratio of the 
inner axisymmetric disk relative to the star is a few per 
cent, so that most of the accreted mass is accumulated 
by the star. The inner disk is dynamically inactive; it 
contributes only to the total gravitational potential and 
its use is necessary to ensure a smooth behaviour of the 
gravity force down to the stellar surface. The other 10% 
of the accreted gas is assumed to be carried away with 
protostellar jets. The latter are triggered only after the 
formation of the central star. 

The basic equations of mass and momentum transport 
in the thin-disk approximation are 

8T 

= -V p -(SO, (1) 



where £ is the mass surface density, V — J Z Z Pdz is 
the vertically integrated form of the gas pressure P, 
Z is the radially and azimuthally varying vertical scale 
height, v p = v r r + v<p4> is the velocity in the disk plane, 
9p = 9rf" + 9it>4> is t ne gravitational acceleration in the 

disk plane, and V p = rd/dr + <f>r~ x d j d<\> is the gradient 
along the planar coordinates of the disk. The gravita- 
tional acceleration g p includes the gravity of a central 
point object (when formed), the gravity of the inner in- 
active disk (r < 5 AU), and the self- gravity of a circum- 
stellar disk and envelope. The latter component is found 
by solving the Poisson integral using the convolution the- 
orem. The disk is pierced with a magnetic field, which 
has only the vertical component B z in the disk and both 
the vertical and planar (B p = B r f + B^<p) components 
at the top and bottom surfaces of the disk. Because of 
our assumption of flux freezing and a spatially uniform 
flux-to-mass ratio a m , the vertical and planar magnetic 
field components are easily determined from the relations 
B z = a m 27rG' 1 / 2 S and B p = — a^gp/G 1 ^ 2 , respectively 
(see iVorobvov fc B asu 2006, for details). 
The viscous stress tensor II is expressed as 

n = 2£i^Vu-i(V-u)e^ , (3) 

where Vi> is a symmetrized velocity gradient tensor, e 
is the unit tensor, and v is the kinematic viscosity. The 
compon ents of (V-II)^ in polar co ordinates (r, cf>) can be 
found in IVorobvov &: Basul (|2009al ). We make no specific 
assumptions about the source of viscosity and parame- 
terize the magnitude of kinematic viscosity using a usual 
form of the a-prescription 

v = a c s . cS Z, (4) 

where c 2 cff = jV/T, is the square of the effective sound 
speed and the ratio of specific heats is set to 7=1.4. We 
note that c Sje ff in this formula is a dynamically varying 
quantity and is calculated at every time step during the 
evolution. 

In this paper, we use a spatially and temporally uni- 
fo rm a = 0.005. This choice is based on the recent work 
of IVorobvov fc Basul (|2009a[ ). who studied numerically 
the secular evolution of viscous and self-gravitating disks. 
They found that if circumstellar disks around solar-mass 
protostars could generate and sustain turbulence then 
the temporally and spatially averaged a should lie in the 
10~ 3 — 10~ 2 limits. Smaller values of a (< 10~ 4 ) have 
little effect on the resultant mass accretion history, while 
larger values (a > 10 _1 ) destroy circumstellar disks dur- 
ing less than 1.0 Myr of evolution and are thus inconsis- 
tent with mean disk lifetimes of the order of 2-3 Myr. 

Equations ([]]) and © are closed with a barotropic 
equation that makes a smooth transition from isother- 
mal to adiabatic evolution at £ = E cr : 

P = C 2 S + c s 2 S cr ^y, (5) 

We note that in this equation (and further in the text) 
c 2 = RTq/[i is the sound speed at the beginning of nu- 
merical simulations, which depends on the adopted initial 
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gas temperature To of a cloud core and the mean molec- 
ular weight p. The value of E cr is calculated during the 
numerical simulations as E cr = m^p n CT 2 Z, where the 
critical volum e number density n cr is set to 10 11 cm~ 3 
(|Larsonll2003| ) and the mean molecular weight is set to 
2.33. The scale height Z is calculated using th e assump- 
tion o f vertical hydrostatic equilibrium (see IVorobvovl 
l2009al) . We note that Z is an increasing function of ra- 
dius, which makes E cr to increase with radius as well. In 
practice, this means that the inner disk regions are signif- 
icantly warmer than the outer regions, since the optically 
thick regime in the inner regions is achieved at lower E. 
This in turn impedes the development of gravitational 
instability and fragmentation in the inner disk, in agree- 
ment with more sophisticated numerical simulations and 
theoretical predictions th at directly solve fo r the energy 
balance equation (see e.g. lBolev et alJfeOOijf ). 

Equations (Q]) and ([2]) are solved in polar coordinates 
(r, <fi) on a numerical grid with 256 x 256 grid zones. We 
use the method of finite differences with a time-explicit 
solution procedure similar in methodology to the ZEUS 
code, with the advection and the force terms treated sep- 
arately using the operator-split method. Advection is 
performed using the second-order van Leer scheme. A 
small amount of artificial viscosity that spreads shocks 
over two grid zones is added. However, the resulted ar- 
tificial viscosity torques are negligible in comparison to 
gravitational and a- viscosity ones. The radial points are 
logarithmically spaced. The innermost grid point is lo- 
cated at r sc = 5 AU, and the size of the first adjacent 
cell varies in the 0.12-0.17 AU range depending on the 
cloud core size. 

3. INITIAL CONDITIONS 

We start our numerical simulations from starless, grav- 
itationally bound cloud cores, which have surface den- 
sities E and angular velocities fi typical for a co llaps- 
ing, a xisymmetric, magnetically supercritical core (jBasul 
1997) 




where f2o is the central angular velocity, ro is the ra- 
dius of central near-constant-density plateau defined as 
r = vC4cs/(7rGE ). We note that when r 3> ro, the 
resultant gas surface density becomes E = s/Ac 2 / (nGr) . 
This profile can be derived from the following gas vol- 
ume density distribution p — Ac 2 / (2irGr 2 ) , if the latter 
is integrated in the vertical direction assuming a local 
vertical hydrostatic equilibrium, i.e. p = E/(2Z) and 
Z = Cg/(7rGE). Therefore, our initial density configu- 
ration is characterized by a factor of A positive density 
enhancement as compared to th at of the s ingular isother- 
mal sphere psis = 

c 2 J{2nGr 2 ) (jShulll977ft . 
Cloud core are characterized by the ratio of rotational 
to gravitational energy /? = -Erot/I-Egravl, where the rota- 
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tional and gravitational energies are calculated as 

^out 1'out 

E rot — 2tt J ra c T,rdr, E glav — — 2ir J rg r T,rdr. 

(8) 

Here, a c = fi 2 r is the centrifugal acceleration, and r out 
is the outer cloud core radius. Cloud cores are initially 
isothermal, with the gas temperature varying between 
To = 10 K and 18 K, depending on the model. 

We present results from 10 sets of models, the param- 
eters of which are summarized in Table [TJ Every model 
set has a distinct ratio (3, the adopted value s of w hich 
lie within the limits inferred by ICaselli et al.l (|2002T) for 
dense molecular cloud cores, /? = (10 -4 — .07). We 
note t hat by construction j3 = 0.91^o r o/ c s (|Vorobvovl 
2009b). In addition, each set of models is characterized 
by a distinct ratio r out /ro in order to generate gravita- 
tionally unstable truncated cores of similar form. For 
most model sets, we choose r ou t/ro = 6.0. 

Every model set contains several individual models 3 
that are characterized by distinct masses M c \, outer radii 
r out , and central angular velocities Qq, but have equal 
/3 and r out /ro. The actual procedure for generating a 
specific cloud core with a given value of f3 is as follows. 
First, we choose the outer cloud core radius r out and find 
ro from the condition r out /ro = const. Then, we find 
the central surface density Eo from the relation ro = 
\f~Ac 2 1 (-/tGEo) and determine the resulting cloud core 
mass M c \ from Equation ([6]) . Finally, the central angular 
velocity ^o is found from the condition (3 = ^lfl^r 2 , / c 2 . 

In total, we have simulated numerically the time evolu- 
tion of 106 cloud cores spanning a range of initial masses 
between 0.06 Mq and 3.9 Mq. We have split this mass 
interval into 14 bins of equal size dM c \. The resulting ini- 
tial cloud core mass function dN/dM c \ = M c 7 m is shown 
in Figure Q] by filled circles. The dashed lines present the 
least-squares best fits to low-mass (0.06 M < M c \ < 
0.6 Mq, left line) cloud cores and intermediate- and 
upper-mass (0.6 Mq < M c \ < 3.9 M s , right line) cloud 
cores. The corresponding exponents are m = 0.5 ±0.1 
and m = 1.6 ± 0.2, respectively. Our model cloud core 
function is similar t o that inferred fr om nearby star- 
forming regions (e.g. lEnoch et al.ll2006T ). The masses of 
our model cloud cores have been selected from a mass 
distribution that mimics the observed core mass distri- 
bution in order to perform averages of the model pa- 
rameters (such as disk masses, Class 0/1 lifetimes, etc.) 
weighted with the cloud core mass function. 

4. CLASSIFICATION SCHEMES 

Modern classification schemes of young stellar objects 
are designed to distinguish between the four main physi- 
cal phases of their evolution: the cloud core collapse and 
formation of a protostar and a disk (Class and Class I), 
envelope clearing and disk accretion onto the star (Class 
II), and disk dissipation and planet formation (Class III). 
The spectral energy distribution is often used to re late 
a YSOs to a particular class (see lEvans et al.l 120091 for 
a thorough review). For instance. lLadal ()1987f ) used the 
spectral index n (defined by vF v cx v n , where F v was 

3 The exact number of these models is specified in the last col- 
umn of Tabic [1] 
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TABLE 1 

Model parameters 



Set 


P 


n 




ro 


Af c i 


r ou t Ad 


To 


A 


Q-m 


N 


1 


2.17 x 10 -3 


0.57 - 3.30 


583 


- 3360 


0.34 - 1.9 


6 


10 


2 





8 


2 


3.20 x 10~ 3 


0.57 - 5.60 


ill 


-4110 


0.24 - 2.4 


6 


10 


2 





12 


2T 


3.20 x 10" 3 


1.30 - 10.5 


111 


- 3360 


0.43 - 3.5 


6 


18 


2 





11 


2A 


3.20 x 10" 3 


0.98 - 9.6 


342 


-3360 


0.40 - 3.9 


6 


10 


8 





12 


2MF 


3.20 x 10" 3 


0.57 - 5.7 


411 


-4110 


0.24 - 2.4 


6 


10 


2 


0.3 


12 


2E 


3.20 x 10~ 3 


1.0 - 5.7 


411 


- 2230 


0.5 - 2.8 


12 


10 


2 





9 


3 


6.90 x 10" 3 


1.11 - 16.6 


206 


- 3086 


0.12 - 1.8 


6 


10 


2 





12 


4 


9.90 x 10" 3 


1.10 - 25.0 


164 


- 3771 


0.095 - 2.2 


6 


10 


2 





10 


5 


1.34 x 10~ 2 


1.33 - 35.0 


137 


- 3600 


0.079 - 2.1 


6 


10 


2 





10 


6 


2.22 x 10~ 2 


1.56 - 45.0 


102 


- 3940 


0.06 - 2.3 


6 


10 


2 





10 


inces are 


in AU, angular velocities in 


km s 


_1 pc~ 1 


masses in M 


3, initial 


gas temperatures Tq 



m=0.5 



\» • 

m=1.6 \ 

\ 



0.1 0.2 0.5 1 2 

Cloud core mass, M d (M@) 

Fig. 1. — Initial cloud core mass function dN/dM c i = M^ m for 
106 model cloud cores organized in 14 bins of equal size dM a \. The 
parameters of the cloud cores are listed in Table [T] The dashed 
lines present the least-squares best fits to low-mass cloud cores 
(0.06 Mq < M c i < 0.6 Mq, left line) and intermediate- and upper- 
mass cloud cores (0.6 Mq < M c i < 3.9 Mq, right line). 



the infrared flux density at frequency u) to distinguish 
between Cla s s I, C lass II, and Class III objects. Later, 
lAndre et al.l (|1993f) introduced the ratio of submillime- 
ter to bolometric luminosity £ S ubmm/-kboi to split the 
early embedded stage into the Class and Class I phases 
based on the mass reservoir remaining in the envelope. 
This separation was originally motivated by the discov- 
ery of very young protostars in collapsing cloud cores 
that were previously being thought of as starless. But it 
also makes physical sense because the Class phase can 
be identified with a period of temporally increasing Lboi 
and Class I phase with a later per iod of decreasing Lbo i 
(IVorobvov fc Basull2005l) . Finally. iMvers fc Laddl (Il993ft 
suggested the use of a bolometric temperature Tboi, de- 
fined as the temperature of a blackbody with the same 
flux-weighted mean frequency as the actual spectral en- 
ergy distribution. 

In our case, it is difficult to use any of the aforemen- 
tioned schemes due to a simplified treatment of the ther- 
mal physics in our numerical simulations. Moreover, the 
large inclination dependence introduced by the outflow 
cavities prevents both Tboi and Lboi/^submm from be - 
ing good evolutionary indicators ([Dunham et al.ll2010D . 
Therefo re, we use the classi fication breakdown suggested 
also by lAndre et al.l (|1993| ) and based on the mass re- 
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Fig. 2. — Azimuthally averaged profiles of the gas surface density 
(top row), radial velocity (middle row), and azimuthal velocity 
(bottom row) for three models from model set 2. The initial cloud 
core masses are indicated in the top row. Four different evolution 
times (since the formation of the protostar) are shown with solid 
(0.1 Myr), dashed (0.2 Myr), dash-dotted (0.3 Myr), and dash- 
dot-dotted (1.0 Myr) lines. The dotted horizontal lines in the top 
row mark the two adopted values of the critical gas surface density 
S d2e for the disk to envelope transition. The dotted line in the 
bottom row presents the angular velocity in the gravitational field 
of a 4.0 Mq star. 

maining in the envelope 



Class M cnv > 0.5M c i 

Class I 0.1M c i < M cm < 0.5M c i 

Class II M cnv < 0.1M cl . 



(9) 



According to this scheme (hereafter, AWTB scheme), 
transition between Class and Class I objects occurs 
when the envelope mass M onv decreases to half of the 
initial cloud core mass M c \. The Class II phase ensues 
by the time when the infalling envelope nearly clears and 
its total mass drops below 10% of the initial cloud core 
mass. Of course, we acknowledge that these numbers are 
somewhat arbitrary and the use of other classification di- 
agnostics may introduce a systematic shift to our derived 
lifetimes, e spectral energy distribution. 

In order to use the AWTB scheme, we need to know 
which part of our numerical grid is occupied by the disk 
and which part belongs to the infalling envelope at any 
time instance of the evolution (note that we know accu- 
rately the stellar mass). This turned out to be a non- 
trivial task. An intuitive step of calculating the rota- 
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tion profile and separating between Keplerian and sub- 
Keplerian regions of the grid sometimes failed to work 
due to large radial motions occasionally present in the 
disk. Besides, the disk is not exactly Keplerian due to a 
noticeable contribution from its self-gravity. Therefore, 
we have developed a hybrid method that makes use of the 
radial velocity field in the envelope and a critical surface 
density for the disk-to-envelope transition. 

Figure[2]illustrates our idea and shows the azimuthally- 
averaged gas surface density E (top row), radial veloc- 
ity v r (middle row), and azimuthal velocity (bottom 
row) distributions as a function of radial distance r for 
three representative models from model set 2. The ini- 
tial cloud core masses M c \ are shown in the top row to 
distinguish between the models. For each model, the ra- 
dial profiles at four different evolutionary times (elapsed 
since the protostar formation) are plotted with the solid 
(0.1 Myr), dashed (0.2 Myr), dash-dotted (0.3 Myr), and 
dash-dot-dotted (1.0 Myr) lines. From the radial distri- 
bution of E it is seen that the disk is compact and dense 
upon formation but later it spreads radially outward due 
to the action of viscou s torques, which are posi tive in 
the outer disk regions (|Vorobvov fc Basul 12009a) ) . The 
disk has a near-Keplerian rotation profile as can be seen 
from the comparison of TJ^ with the dotted line, the latter 
showing the angular velocity in the gravitational field of a 
4.0 M & star. It is worth noting that the evolution of the 
unmagnetized models seems to approach a self-similar 
behavior in the infalling envelope. This is evident in the 
M = 1.0 Mq model, which radial velocity and surface 
density distributions approach a free-fall profile propor- 
tional to r -1 / 2 . 

In our numerical simulations, the disk occupies the in- 
ner regions (small r), while the envelope — the outer re- 
gion of the computational grid (large r) . A large drop in 
v T (by absolute value) correspond to the radial position 
where the infalling envelope lands onto the disk (the so- 
called accretion shock). The shock position propagates 
radially outward, reflecting the disk expansion. It is also 
evident that the shock front becomes less sharp with in- 
creasing cloud core mass, reflecting an increasingly com- 
plex spiral structure of massive disks (recall that Figure^] 
presents azimuthally averaged profiles). It is seen that, 
in general, the envelope is characterized by large neg- 
ative v r due to gravitationally driven collapse and low 
E due to gas depletion. On the other hand, the disk 
has much smaller v r due to near-(but never complete)- 
centrifugal balance and larger E due to mass accumula- 
tion. The outer parts of the disk expand and the overall 
disk density decreases with time due to ongoing angular 
momentum redistribution. In general, we find that mas- 
sive disks are not in steady state but exhibit substantial 
radial pulsations due to the time-dependent strength of 
gravitational instability in the disk. The gas surface den- 
sity profile declines with radius but some nonmonotonic 
behavior is evident for upper-mass models due to frag- 
mentation in the outer parts of the disk. 

Having analyzed these data, we adopt the following 
procedure to distinguish between the infalling envelope 
and the disk. We march from the outer computational 
boundary toward the inner one (in the direction of de- 
creasing r) and determine the radial distance at which 
the following two criteria are fulfilled simultaneously: 
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Fig. 3. — Time evolution of the stellar mass (solid lines), envelope 
mass (dashed lines) and disk mass (dotted and dash-dotted lines) 
for four models from model set 2. The initial cloud core masses M c i 
are indicated in each panel. Two disk masses based on the value of 
the critical gas surface density £d2c are shown by the dash-dotted 
lines (Sd2c = 0-1 g cm -2 ) and dotted lines E^e = 0.01 g cm -2 . 
The vertical dotted lines mark the onset of the Class I phase (left 
line) and Class II phase (right line). 

E > Ed2o and v r > 0. This radial distance is set to 
represent the disk's tentative outer radius r<j (at a given 
time instance of the evolution). For the critical transi- 
tional density Ed2e, we choose 0.01 g cm -2 or 0.1 g cm" 2 
shown in Figurc[2]with horizontal dotted lines. All mate- 
rial that is located within r < is considered to belong 
to the disk. This procedure is repeated each time we 
want to determine the disk and envelope masses and it 
has proven to be most accurate in our circumstances. 

Figure [3] presents the integrated disk (dotted and dash- 
dotted lines) , envelope (dashed lines) , and stellar masses 
(solid lines) as a function of time elapsed since the begin- 
ning of numerical simulations in four representative mod- 
els from model set 2. The initial cloud core masses M c \ 
are indicated in each panel. Two different disk masses are 
shown based on the critical density Ed2e = 0.01 g cm -2 
(dotted lines) and Ed2c = 0.1 g cm~ 2 (dash-dotted lines). 
The vertical dotted lines mark the onset of the Class I 
phase (left line) and Class II phase (right line). It is seen 
that the appearance of the star occurs at progressively 
later evolution times with increasing cloud core mass. 
This is due to the fact that the free-fall time Tg of a 
cloud core is linearly proportional to its mass (see the 
Appendix). Indeed, in the M c i=0.38 Mq model, the star 
appears at t = 0.026 Myr, whereas in the M c \—1.9 Mq 
model the star appears at t = 0.129 Myr, which is in 
good agreement with what is expected from the linear 
proportionality between and M c \. It is also seen that 
the difference between disk masses based on the two val- 
ues of Ed2c is mostly indistinguishable, apart from sev- 
eral episodes when some noticeable mismatch is evident. 
These episodes occur only in the early evolution (and 
in massive disks) and are associated with a transitory 
disk expansion after the accretion of massive fragments 
(for med due to disk fragment ation) by the central proto- 
star (|Vorobvov k. Basull2006[) . 

Figure [3] also shows that the disk mass is always 
lower than that of the protostar. This is typical of 
self-gravitating d isks in which sel f -gravity is compute d 
self-consistently (|Vorobvovl l2009at iKratter et alj|2009l) . 
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rather than being calculated using an effective vis- 
cosity proport i onal t o the Toomre Q-parameter as in 
iLin fc Pringlel (|1990() . In the latter case, the stellar 
mass may be substantially underestimated in models 
with large disk-to-star mass ratios, whereas the disk mass 
agrees better with that obtain ed from self-consi stent cal- 
culations of disk self-gravity (|Vorobvovl l2009cf) . A net 
result of this discrepancy is that the disk mass starts to 
considerably exceed that of the star soon after its for- 
mation, which contradicts observations. We note that 
the time difference between the star and disk formation 
events is somewhat overestimated in our numerical sim- 
ulations due to the use of the sink cell. The rate at 
which the envelope mass declines with time is clearly de- 
pendent on M c \ — the envelope depletes much slower in 
models with greater M c \. 

5. LIFETIMES OF THE EMBEDDED PHASE 

In this section, we calculate lifetimes of the Class 
and Class I phases based on the AWTB criterion (see 
Equation [§])• Depending on the adopted value of the 
critical density Ed2e for the disk-to-envelope transition 
(0.01 g cm~ 2 or 0.1 g cm -2 ), the resulting lifetimes 
may differ by a factor of unity. In order to minimize 
the uncertainty, we use the arithmetically averaged val- 
ues. For model set 6, however, we use lifetimes based 
on Ed2e = 0.01 g cm -2 . This is because model set 6 
produces most massive and extended disks (due to con- 
siderable rotation rates of the corresponding cloud cores) , 
which are best described by the smaller value of Ed2c- 

For the adopted initial (isothermal) gas surface den- 
sity distribution of the form £ oc r , one may expect 
that Class and Class I lifetimes (tqo and tci, respec- 
tively) scale linearly with cloud core mass M c \, Indeed, 
for this case, the enclosed mass is proportional to radius 
and therefore the crossing-time t CT — r/c s of an infalling 
gas shell is roughly proportional to the enclosed mass. 
A more rigorous derivation given in the Appendix shows 
that the free-fall time Tg is linearly proportional to M c \, 
inverse proportional to the initial density enhancement 
amplitude A, and scales with the initial gas temperature 
as T~ 3 / 2 . If the ratios tqo/ts and tci/t^ are indepen- 
dent of cloud core mass (which is not always the case, 
see Figure [7]), then a linear relation between lifetimes 
and cloud core masses may also be anticipated. 

However, the Class 0/1 lifetimes are often inferred from 
observations using stellar characteristics, and not those 
of cloud cores. Therefore, it is desirable to relate life- 
times with observationally meaningful stellar properties 
such as stellar masses. The existence of a linear depen- 
dence between lifetimes and stellar masses, similar to 
that anticipated for cloud core masses, is less obvious a 
priori. First, the stellar mass varies significantly dur- 
ing the Class phase and, to a lesser extent, during the 
Class I phase (see Figure [3]). It could be observation- 
ally difficult to know the exact evolutionary instance of 
the forming star, i.e., whether the star is near the be- 
ginning of, say, the Class phase or in the midway of 
it. Second, not all infalling mass goes into the star. A 
significant fraction lands onto the disk, and the disk-to- 
star mass ratio also varies si gnificantly during the early 
evolution (Figure [31 see also IVorobvovl (|2009aD ). These 
facts suggest that a linear correlation between lifetimes 
and stellar masses can be theoretically anticipated only 



if the latter are the final masses. On the other hand, 
using the final stellar masses in the lifetime-stellar mass 
dependence may be misleading, since observations would 
provide a whole spectrum of masses at various evolution 
instances of the Class and Class I phases, and not only 
at the end of them. 

To complicate things even further, the lifetime-stellar 
mass correlation may be influenced by other cloud core 
parameters such as magnetic fields and the initial rate 
of rotation, the effects of which are difficult to eval- 
uate theoretically. And finally, although the adopted 
S oc r _1 scaling of the initial gas s urface density with ra- 
dius is borne out o bservationally (iBacmann et al.l [20001 : 
iDapp fc Basu| [2009) and theoretically (|Basulll997lh there 
is no guarantee that other scaling laws could not be re- 
alized in nature. This justifies a detailed numerical in- 
vestigation into to the problem. 

5.1. Lifetimes versus stellar masses 

In this section, we present numerically derived rela- 
tions between Class and Class I lifetimes and corre- 
sponding stellar masses, M^co and M^ci- Since stel- 
lar masses quickly increase with time (see Figure [3]), we 
choose to consider five typical stellar masses based on the 
ratio M cnv /M c \. This ratio is useful because it gradually 
decreases with time and can be used as a tracer of the 
early evolution of a protostar. In particular, we consider 
three stellar masses derived when the ratio M env /M c \ 
becomes equal 0.9, 0.7, and 0.5. These stellar masses 
are meant to represent typical ones in the beginning, in 
the midway, and at the end of the Class phase, re- 
spectively. For the typical stellar masses in the Class I 
phase, we choose three values based on M env /M c i=0.5, 
0.3, and 0.1. We note that we use the same stellar masses 
at the end of the Class phase and at the beginning of 
the Class I phase, but the corresponding lifetimes are 
different. Of course, these five stellar masses cannot rep- 
resent the whole possible spectrum of stellar masses that 
can be observationally detected in the embedded phase 
of star formation. Nevertheless, they can help to reduce 
a possible bias towards a particular evolution stage and 
can help to make our lifetime-stellar mass relation more 
observationally meaningful. 

Figure H] presents our model lifetimes (in Myr) of the 
Class phase (top panels) and Class I phase (bottom 
panels) versus stellar masses (in M Q ) for six model sets. 
More specifically, model set 1 is plotted by red triangles- 
up, model set 2 — by blue circles, model set 3 — by green 
squares, model set 4 — by plus signs, model set 5 — by 
pink diamonds, and model set 6 — by brown triangles- 
down. Stellar masses in each panel are derived at specific 
evolutionary times as specified by the ratio of envelope 
mass to cloud core mass M env /M c i. In particular, the 
upper-left panel employs stellar masses obtained at the 
beginning of the Class phase, respectively The upper- 
middle and upper-right panels use stellar masses derived 
in the midway and at the end of the Class phase. The 
same age segregation is used in the bottom row of panels 
only now the stellar masses are related to the Class I 
phase. 

Each symbol of same color and shape within a given 
set of models represents an individual object, which has 
formed from a cloud core of distinct mass, rotation rate, 
and outer radius. 
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Fig. 4. — Lifetimes of the Class phase (top row) and Class I phase (bottom row) versus stellar masses calculated at the beginning (right 
column), in the midway (middle column), and at the end (right column) of each evolution phase. The data for six model sets (see Table [TJ 
are shown. The solid lines present the least-squares fits to these data. The ratio M eIlv /M c \ of the envelope mass to the initial cloud core 
mass and the exponent m of the least-squares fit are indicated in every panel. 



There are two interesting features in Figure 2] that de- 
serve particular attention. First, and most important 
one, is a notable dependence of the derived lifetimes on 
the mass of the central object. The least-squares best 
fits (solid lines) yield a near-linear correlations between 
the Class lifetimes and stellar masses and a somewhat 
steeper correlation for the case of the Class I lifetimes. 
The corresponding exponents m are indicated in every 
panel. Second interesting feature in Figure|4]is seen along 
the line of constant M* — there is a trend for models with 
greater ratio of rotational to gravitational energy /3 to 
have greater Class 0/1 lifetimes. This effect is especially 
prominent for the Class I objects and is not unexpected. 
Cloud cores that have more energy stored in rotation are 
expected to retain the envelope for a longer time due to 
an increased centrifugal force, thus prolonging the em- 
bedded phase of star formation. 

The correlations shown in each panel of Figure [4] may 
not be observationally meaningful when taken separately, 
because they are derived for stars of specific age. As 
noted above, observations would not be biased towards 
any specific instance of protostellar evolution but arc 
likely to detect a whole spectrum of stellar ages. To 
make our predictions more observationally significant, we 
plot in Figure [5] the Class (top) and Class I (bottom) 
lifetimes versus all stellar masses, namely, those at the 
beginning, in the midway, and at the end of each evo- 
lution phase. In other words, in the top/bottom panel 
of Figure [5] we combine the three top/bottom panels of 
Figure |H The least-squares best fits (solid lines) yield 
the following relations between Class (too) and Class I 
(tci) lifetimes (in Myr), from one hand, and the corre- 



sponding stellar masses (in Mq), from the other hand 



rco = 0.18 M, 
T C i = 0.30 M. 



0.8±0.05 

*,C0 ' 

1.2±0.05 

*,CI 



(10) 
(11) 

It is evident that a somewhat steeper than linear cor- 
relation between tqi and M^ci persists irrespective of 
the stellar age, i.e., independent of whether we measure 
stellar masses in the beginning, in the midway, or at the 
end of the Class I phase. This is a direct consequence 
of the fact that a substantial fraction of the final stellar 
mass has already been accumulated by the beginning of 
the Class I phase and the subsequent increase in M^ci 
is not significant (see Figure [3|). On the other hand, the 
Class lifetimes demonstrate a noticeably weaker cor- 
relation with stellar masses. This is because the latter 
vary substantially during the Class phase, which effec- 
tively spreads stellar masses along the horizontal axis in 
Figure [3 thus lowering the exponent in the least-squares 
fit. 

This correlation between lifetimes and stellar masses 
suggests that the disagreement between the observation- 
ally inferred Class and Class I lifetimes of YSOs may 
in part be due to variations in the initial mass function 
or due to observational biases toward a particular mass 
band. To illustrate this idea, we calculate mean lifetimes 
for objects in our entire mass band (0.008-0.85 M Q and 
0.03-1.24 M Q for Class and Class I objects, respec- 
tively) and for objects in the truncated mass band (0.2- 
0.85 M Q and 0.2-1.24 M© for Class and Class I objects, 
respectively), with the latter case being confined to inter- 
mediate and upper-mass stars only. The resultant mean 
lifetimes of the Class objects are (too) = 0.044 Myr 
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Fig. 5. — Lifetimes of the Class phase (top) and Class I phase 
(bottom) versus stellar masses, the latter including those calculated 
at the beginning, in the midway, and at the end of each evolution 
phase. The top /bottom panels of the figure are obtained by com- 
bining the top/bottom rows of Figure [4] Solid lines present the 
least-squares best fits with exponents m. 

and (r co )tr = 0.086 Myr, respectively. In the case of 
the Class I objects, we obtain (rci) = 0.09 Myr and 
(Tci)tr =0.15 Myr. It is evident that the neglect of the 
low-mass objects causes almost a factor of two overesti- 
mate of the Class 0/1 lifetimes. In spite of the fact that 
Too show a somewhat weaker correlation with M^co, 
the effect of truncation is even greater than in the case 
of Class I lifetimes due to systematically smaller stellar 
masses in the Class phase. 

5.2. The effect of varying initial conditions and 
magnetic fields 

In the previous section, we have considered model cores 
that differ in the adopted values of rotational to grav- 
itational energy /3, initial mass M c \, and radius r out . 
However, cloud cores may vary in other aspects such as 
initial gas temperatures, truncation radii, initial density 
enhancements, and strengths of magnetic fields. In ad- 
dition, the shape of the initial gas surface density profile 
may differ from the assumed scaling K oc r' 1 . It is com- 
putationally prohibitive to consider the effect of vary- 
ing initial conditions on every model core of Section 15.11 
Therefore, we have chosen model set 2 as the reference set 
to estimate the effect that the varying initial conditions 
may have on the lifetime of the embedded phase. 

Figure IS] shows our model lifetimes of the Class phase 
(top panels) and Class I phase (bottom panels) versus 
stellar masses, the latter derived as described in the pre- 
vious section, for five model sets. In particular, red cir- 
cles correspond to model set 2, blue squares correspond 



to model set 2T, which is characterized by a higher ini- 
tial gas temperature of the cloud cores T = 18 K. Green 
triangles-up present model set 2E, which is character- 
ized by extended cloud cores with r ou t/^o = 12 (r out is 
increased). Pink diamonds show model set 2A, in which 
cloud cores are characterized by the initial density en- 
hancement amplitude A = 8, and yellow triangles-down 
correspond to model set 2MF with the flux-to-mass ra- 
tio a m — 0.3. Finally, plus signs present model set 2C, 
which comprises models with spatially constant surface 
density and angular velocity distributions. The latter are 
obtained using the parameters of the reference model set 
and setting E = E(r = r out ) and Q = ft(r = r out ) in 
Equations ^ and (0, respectively. In order words, sur- 
face densities and angular velocities in model set 2C are 
set equal to those at the boundary in the corresponding 
models of the reference set. The ratio /3 = 2.3 x 10~ 3 
is the same for all models in this section in order to ex- 
clude a mild dependence of the lifetimes on (3 discussed 
in Section [5TT1 

It is evident that the varying initial conditions may 
have a profound effect on the duration of the embed- 
ded phase. In particular, higher initial gas temperatures 
and density enhancement amplitudes reduce the result- 
ing Class 0/1 lifetimes by factors of 2.5-5.5 for objects 
with equal mass M*. This is not surprising, since both 
effects are expected to increase the mass accretion rate 
from the envelope M cnv , thus shortening Class 0/1 life- 
times. In the standard model of inside-out collapse, M env 
is proportional to the cube of the sou nd speed c f and to 
the density enhancement amplitude A (|Shulll977t ). In our 
case, the increase in gas temperature is AT =1.8 and 
Mcnv is expected to increase by a factor of 1.8 3 / 2 = 2.4. 
Indeed, top panels in Figure [6] indicate that the corre- 
sponding Class and Class I lifetimes (blue squares) have 
decreased on average by factors of 2.5 and 3.2, respec- 
tively, as compared to the reference set of models (red 
circles). A somewhat greater drop in tqi than expected 
indicates that the envelope empties in the late embedded 
phase faster than predicted by the standard model. This 
could be due to the influence of the cloud core's finite 
outer boundary and an inward-propagating rarefaction 
wave (cloud cores in the standard model are infinite in 
size). This effect becomes even more pronounced in the 
case of a greater initial density enhancement A — 8 (as 
opposed to A = 2 in the reference set). The correspond- 
ing Class and Class I lifetimes, shown by the pink di- 
amonds in Figure [6l have dropped by factors of 2.6 and 
5.5, respectively. 

On the other hand, YSOs formed from cloud cores with 
greater truncation radii r out /rn = 12 (green triangles-up 
in Figure do not seem to have significantly different 
lifetimes as compared to those of the reference model 
set with rout/ro = 6. This may seem counterintuitive. 
Indeed, more extended cores have a greater mass reser- 
voir and are supposed to deposit their material onto a 
star/disk system for a longer time. However, one has to 
keep in mind that they would also form more massive 
stars. As a result, the lifetimes simply shift toward the 
upper-right corner in every panel of Figure [5] 

YSOs formed from cloud cores with frozen-in magnetic 
fields (a m = 0.3) reveal a complicated behavior. On 
one hand, their Class lifetimes, as shown by yellow 
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triangles-down in the top-left panel of Figure El do not 
seem to be significantly affected. On the other hand, 
Class I lifetimes of the intermediate- and upper mass 
stars increase considerably as compared to the reference 
set of model. The reason why these stars are more af- 
fected than their lower-mass counterparts lies in the way 
that the frozen-in magnetic fields operate in cloud cores. 
Let us consider the outermost layer of gas. As it moves 
toward the center, magnetic tension increases and de- 
celerates the collapse, prolonging the embedded phase. 
Magnetic field lines are attached at infinity to the natal 
giant molecular cloud, which evolves on a much slower 
time scale than the collapsing cloud core. This effec- 
tively means that the field lines are motionless at infinity. 
It now follows that magnetic tension would be stronger 
for larger (and hence more massive) cores, as it takes a 
larger distance for the outermost layer to travel toward 
the disk surface. If we recall that more massive cores 
produce more massive stars, a steeper dependence of tqi 
on -M^ci becomes evident. We note that this effect may 
be either magnified if higher values of a m are present in 
cloud cores or moderated for lower values of a m . In ad- 
dition, ambipolar diffusion and magnetic breaking, not 
included in the present numerical simulations, are ex- 
pected to decrease the effect. 

Finally, the initial distribution of gas surface density 
and angular velocity in cloud cores may also affect the re- 
sulting lifetimes. For the case of cloud cores with initially 
constant E and f2 (plus signs in Figure [6]), the resulting 
lifetimes are greater than those in the reference model 
set, the latter having both S and inverse proportional 
to radius at large radii. The increase is particularly large 



in the case of Class I lifetimes and is caused by a steep 
radial gas surface density distribution that develops in 
the late phases of cloud core collapse. Although cloud 
cores in model set 2C start with constant E, they quickly 
evolve into a configuration that is characterized by a gas 
surface density profile that declines with radius steeper 
than r -1 at large radii. This acts to lower the mass in- 
fall rate from the envelope M = — 27rri> r E as compared 
to that in the reference model, the latter being charac- 
terized by E cx r _1 at r >> rn. The effect is particularly 
strong in the late evolution when gas layers initially lo- 
cated at large radial distances begin to fall in onto the 
disk. 

How do varying initial conditions affect the correlation 
between Class /I lifetimes and masses of the central ob- 
ject discussed in the previous section? To quantitatively 
answer this question would require to run a prohibitively 
large number of models for every model cloud core with 
all possible initial conditions. However, based on the be- 
haviour of model cores considered in this section, we can 
assume that the correlations given by Equations (fT0|) and 
(jllj) arc unlikely to break completely. Indeed, in most 
model sets with varying initial conditions (apart from 
that with non-zero frozen-in magnetic fields) the slope of 
the correlation line is similar to that of model set 2. It 
is thus the spread in the lifetimes that can affect some- 
what the resulting correlation. Just to estimated the 
effect, we have combined Class and Class I data from 
all models into two separate t-M„ plots for each phase 
and found that the resulting least-squares exponents in 
Equations (fTU| and (fTTj) decrease by as much as 0.3. 
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Fig. 7. — Class (black bars) and Class I (grey bars) lifetimes 
versus initial cloud core masses for six model sets as specified in 
each panel. Solid lines show the free-fall time rg out of the outer- 
most gas layer as derived from Equation (I A3 ft ■ 

5.3. Class I lifetime versus Class lifetime 

It is often useful to know how Class and Class I life- 
times are related to each other for a particular YSO. It 
is difficult to extract such information directly from Fig- 
ures |3HHl and Equations (fTCTj) and (ITT1) are not useful be- 
cause they provide data averaged over all model cores. In 
order to analyse the relative lifetimes, we have to choose 
some characteristic of a YSO in order to present tqo and 
rci as a function of this characteristic. We use the initial 
cloud core masses with a practical purpose to verify if the 
lifetimes indeed scale linearly with M c \, as was suggested 
on theoretical grounds at the beginning of Section [3] 

Figure [7] presents the Class (black bars) and Class I 
(grey bars) lifetimes as a function of the initial cloud 
core mass M c \ for model cores from six model sets as 
specified in each panel. Every pair of bars represents 
an individual YSO formed from a cloud core with mass 
M c \. It is evident that rci are systematically greater 
than rco, except for model set 2 A for which rci < tqo- 
We note that the latter case is characterized by a rather 
large initial density enhancement A — 8 and need specific 
conditions (such as strong shock waves) to be realized. It 
is also seen that the ratio tci/tco varies little with M c i, 
with a notable exception for models with magnetic fields 
(model set 2MF) and for models with large rotation rates 
(model set 5 and 6). We have already discussed above 
the case with magnetic fields. As to the models with 
greater values of the rotational to gravitational energy /3, 
they have more centrifugal support against gravity, which 
effectively prolongs the Class I phase in the upper-mass 
(and spatially extended) cloud cores. Excluding model 
set 2A and the upper mass models in model sets 5, 6, and 
2MF, the mean ratio of the Class I to the Class lifetime 
is 7ci/tco ~1.5— 2.0. For the aforementioned upper-mass 



models, the corresponding ratio may be as high as 3.0- 
4.0. 

The solid lines in Figure[7]show the free-fall time Tff jOU t 
of the cloud's outermost layer as a function of M c \. The 
free-fall time is linearly proportional to M c \ (see equa- 
tion (|A3[) in the Appendix) . By comparing lifetimes with 
the free-fall time, one can notice that both the Class 
and Class I lifetimes in model sets 2, 2T, and 2A are also 
linearly proportional to M c \, as was anticipated from sim- 
ple theoretical grounds in Section[SJ However, the Class I 
lifetimes of the upper-mass model cores in model sets 5, 
6, and 2MF show a significant deviation from the lin- 
ear scaling typical for the low- and intermediate-mass 
models. This effect was not envisioned from simple the- 
oretical considerations. 

6. ENVELOPE DEPLETION RATES 

In semi-analytic models of cloud core collapse and 
star/disk formation, it is often needed to know 
the rate at which the envelope material is accreted 
onto the star/disk system. We determine this rate 
by calculating the envelope depletion rate M cnv = 
- (M env (t + At) - M env (t)) /At, with the time step At 
set to 10 3 yr. In the earliest stages of evolution before 
the disk formation, M cnv would represent the mass ac- 
cretion rate onto the star, but when the disk forms M env 
would represent the mass accretion rate onto the disk. 

The solid lines in Figure [8] show the envelope depletion 
rate versus time elapsed since the beginning of numeri- 
cal simulations for the same four representative models 
as in Figure [3] The vertical dotted lines mark the onset 
of the Class I (left line) and Class II (right line) phases. 
The time evolution of M env can be split into two distinct 
modes: a shorter period of near-constant depletion rate 
and a longer period of gradual (and terminal) decline of 
M cnv . It is seen that the boundary between these two 
modes lies near the vertical dotted line that separates 
the Class and Class I phases of star formation in ev- 
ery model. A sim il ar be havior of M cnv was reported by 
IVorobvov fc Basul (|2005t ) for the case of the spherically 
symmetric collapse of truncated cloud cores. 

In the first mode, the accretion process is conceptually 
similar t o that desc ribed by the single isothermal sphere 
solution (|Shulll977t) , wherein the material is falling freely 
onto the star/disk system and the accretion rate is con- 
stant and proportional to c^/G. The second mode en- 
sues when a rarefaction wave propagating inward from 
the core's outer boundary hits the disk surface. From 
this moment on, the accretion process is no more self- 
similar, the radial gas density profile starts to steepen 
and the ongoing depletion of gas in the envelope causes 
a gradual (and terminal) decline of M cm . 

We note that M cnv may exhibit some short-time fluctu- 
ations in the second mode due to temporary disk expan- 
sions and contractions caused by disk radial pulsations 4 . 
When expanding, the disk delivers part of its material to 
the envelope, thus temporarily decreasing the envelope 
depletion rate M cnv . The fluctuation amplitudes grow 
along the line of increasing cloud core masses (and, con- 

4 The animation of this process can be viewed at 
www.astro.uwo.ca/~vorobyov (animations: burst mode of accre- 
tion) 
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Fig. 8. — Envelope depletion rates M cnv (or, equivalently, mass accretion rates onto a star/disk system) versus time elapsed since the 
beginning of numerical simulations for four representative models from model set 2. The corresponding cl oud core masses are shown in 
each panel. The solid lines show the model data, while the thick dashed lines are derived using Equation I I14I) , The dashed-dotted and 
dash-dot-dotted lines show M env as derived using the RMA function (1 1 3 P with rg = Tg , out and Tg = iff 2/3, respectively (see the text for 
details). The vertical dotted lines mark the onset of the Class I (left) and Class II (right) phases. 

was constrained by iBontemps et al.l (|1996l) using the in- 
ferred duration of the embedded phase, and Tg is the 
free-fall time. 

First, we check the utility of the BATC function. When 
one undertakes modeling of the gravitational collapse 
of a cloud core, the value of r c is not known a priori. 
Therefore, it is desirable to link r c with the free-fall time 
T ff = (37I-/32G/9) 1 / 2 . We have experimented with various 
values of Tg , depending on the actual value of p, and rec- 
ommend the following relation r c = Tff iOUt /3, where Tg tOUt 
is the free-fall time of the outermost layer derived in the 
Appendix. The time t in Eqs. (|T2"j) and (fT3")) is counted 
from the formation of the protostar (the accretion rate 
is negligible at the earlier time). The resulting mass ac- 
cretion rate Mbatc is shown in Figure by the thick 
dashed lines. It is evident that the BATC function pro- 
vides an adequate fit to our model accretion rates (solid 
lines). This result is robust and holds for models with 
different values of /? and with different initial conditions. 

We now proceed with considering the utility of the 
RMA function. Since our cloud cores are initially non- 
uniform, we have considered several values of Tg, de- 
pending on the actual value of p. The dash-dotted 



sequently, increasing disk masses), indicating that gravi- 
tational instability lies behind these fluctuations. These 
fluctuations, though having a conceptually different na- 
ture, are similar in beha vior to the spasmodic accretio n 
phenomenon reported by iTassis fc Mouschoviasl (|2005f ) . 

We now want to find a simple and convenient way to 
approximate our model envelope depletion rates. We 
co nsider two different fun ctions: the first one proposed 
by IBontemps et alj (|1996() to account for the observed 
outflow energetics in the embedded phase (hereafter, 
the BATC funct ion) and the second one employed by 
iRice et all ((2009) to mimic the mass infall rates onto the 
star/disk system in their one-dimensional global simula- 
tions of cloud core collapse (hereafter, the RMA func- 
tion) 



Mbatc(<) = — x 

Tr. 



(12) 



^ < t < Tg 

-U 1(M ,( / i { IM£ [l + COS Tg < t < 2Tg 

2rg < t. 

(13) 

Here, t c is some characteristic time, the value of which 
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and dash-dot-dotted lines in Figure [8] show the result- 
ing mass accretion rate Mrma for rg = Tff.out an d 

1/2 

m = Tff,2/3 = (37r/(32Gp 2 /3)) , respectively, where 
the latter value is the free-fall time of a layer located at 
r = 2r out /3 (i.e. at a radial distance equal to 2/3 that of 
the cloud core outer radius) and p 2 /3 is the correspond- 
ing gas volume density. It is obvious that, regardless of 
the value of ts, the RMA function fails to provide an 
acceptable fit to our model accretion rates (solid lines). 
To summarize, we suggest using the following equation 
as the mass accretion rate onto a star/disk system 

M„(()^xe- 3t / T «-, (14) 

Tff,out 

where time t should be counted from the formation of a 
central star and Tff, ut (see Equation |A3| ) is the free-fall 
time of the outermost gas layer of a cloud core with mass 
Mci- 

7. CONCLUSIONS 

Using numerical hydrodynamics simulations, we have 
computed the gravitational collapse of a large set of 
cloud cores. We start from a gravitationally unsta- 
ble, pre-stellar phase and terminate the simulations in 
the Class II phase when most of the cloud core has 
been accreted by the forming star/disk system. We 
have considered cloud cores with various initial masses 
(M cl =0.06-3.9 M ), ratios of the rotational to the grav- 
itational energy (/3=(0.2-2.2)x 10 -2 ), initial gas temper- 
atures (To=10-18 K), truncation radii (r out /r =6-12), 
strengths of frozen-in magnetic fields (a m =0-0.3), ini- 
tial gas density enhancements (A=2— 8), and initial radial 
profiles of the gas surface density and angular velocity 
(E, n oc r _1 and S, ft = const). 

We employ a sophisticated method for distinguishing 
between the infalling envelope and the forming disk in 
our numerical simulations and determine the duration of 
the embedded phase of star formation, adopting a clas- 
sification scheme based on the remaining mass in the 
envelope ([Andre et al.l ll993). We also calculate the en- 
velope depletion rate M cnv (or, equivalently, the rate of 
mass accretion onto the star/disk system) and test the 
utility of two empirica l functions for M env p rovided by 
iBontemps et al.l (|1996t) (BATC function) and lRice et~atl 
d2009|) (RMA function). We find the following. 

1. Class (rco) and Class I (rci) lifetimes corre- 
late with the corresponding stellar mass M^co and 
M^ci, however the scaling is more complex than 
can be expected from simple theoretical grounds 
based on the linear correlation between the free- 
fall time and the cloud core mass. In particular, the 
form of the correlation depends on the cloud core 
properties such as the initial ratio of rotational to 
gravitational energy, initial gas temperature, den- 
sity enhancement amplitude, and the initial radial 
distribution of gas surface density and angular ve- 
locity. In addition, frozen-in magnetic fields with 
constant flux-to-mass ratio can substantially influ- 
ence the Class I lifetimes. The correlation is further 
complicated due to the fact that the stellar mass 
varies considerably during the Class phase. This 
makes the rco-M^co correlation dependent on the 
stellar age as well. 



2. When cloud cores with varying rotation rates, 
masses and sizes (but otherwise identical) are con- 
sidered, Class and Class I lifetimes have sub- 
and super-linear correlations with the correspond- 
ing stellar masses, respectively. These correlations 
extend from sub-stellar masses (~ 0.01 Mq) to 
at least solar masses (~ 1.5 M ) and obey the 
following scaling laws r co = 0.18 M°^ 05 and 
rci = 0.3 Mj; 2 f 05 for Class 0/1 stars of all pos- 
sible ages. This makes the observational determi- 
nation of the mean lifetimes sensitive to the form 
of the initial mass function and/or to instrumen- 
tal biases toward a particular mass band. For in- 
stance, our modeling predicts a mean Class life- 
time of (rco) = 0.044 Myr for stars in the 0.008- 
0.85 M mass range and (rco,tr) = 0.086 Myr for 
stars in the 0.02-0.85 M Q range. In the case of the 
Class I objects, we obtain (tci) = 0.09 Myr and 
(rci)tr = 0.15 Myr for stars in the 0.03-1.24 M Q 
and 0.2-1.24 M© mass ranges, respectively. It is ev- 
ident that the neglect of objects at the lower mass 
end results in almost a factor of 2 overestimate of 
the mean lifetimes. 

3. An increase in the initial cloud core temperature 
and density enhancement amplitude tend to lower 
the Class and Class I lifetimes, whereas cloud 
cores with initially constant gas surface density 
and angular velocity distributions (as compared to 
those with S,f2 oc r" 1 ) have longer Class and 
especially Class I lifetimes. In addition, frozen- 
in magnetic fields may increase the Class I life- 
times, particularly for the upper-mass cloud cores, 
thus steepening the corresponding tci-M^ci rela- 
tion. The net effect of varying initial condition is 
to weaken the aforementioned sub-linear (Class 0) 
and super-linear (Class I) correlations between the 
lifetimes and stellar masses by decreasing the cor- 
responding exponents by as much as 0.3. However, 
more accurate modeling with ambipolar diffusion 
and magnetic braking taken into account is needed 
to accurately assess the magnitude of this effect. 
The outer truncation radius of a cloud core has lit- 
tle effect on the resulting lifetimes. 

4. Most cloud cores give birth to YSOs whose Class I 
lifetimes are longer than those of the Class phase 
by roughly a factor of 1.5-2.0. A notable excep- 
tion are YSOs formed from cloud cores a with large 
initial density enhancement. In the latter case, 
the duration of the Class I phase may actually 
be shorter than that of the Class 0. In addition, 
the upper-mass models M c \ > 1 Mq with frozen-in 
magnetic fields and high cloud core rotation rates 
/? > 10~ 2 may have the tqi/tqq ratios as large as 
3.0-4.0. 

5. The time evolution of M cnv reveals two distinct 
modes: a shorter period of near-constant depletion 
rate and a longer period of gradual (and termi- 
nal) decline of M env . The boundary between these 
two modes lies near the end of the Class phase 
and the beginning of Class I phase. In the later 
mode, M cnv may show short-term fluctuations due 
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to episodic disk expansions and contractions. The 
BATC function may provide an adequate fit to our 
model envelope depletion rates if the characteris- 
tic time of decline of M onv is chosen properly. The 
RMA function fails to provide an acceptable fit to 
our model data, irrespective of the free parameters. 

We emphasize that our lifetimes have been derived 
based on the AWTB classification scheme (jAndre et al.1 
1993) and may change by a factor of unity if other 
schemes are used. For instance, if we rc-definc the bound- 
ary between the Class I and Class II phases in the AWTB 
scheme and assume that the Class II phase begins when 
the envelope mass drops below 5% of the initial cloud 
core mass (in contrast to 10% adopted in our paper), the 
resulting Class I lifetimes in model set 2 increase by a 



factor of 1.5. However, the derived trends and correla- 
tions between the lifetime of the embedded phase, from 
one hand, and the stellar masses, from the other hand, 
are expect to stay valid irrespective of the classification 
scheme used. 
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APPENDIX 

RELATION BETWEEN THE CLOUD CORE MASS AND THE FREE-FALL TIME 

If the initial gas volume density distribution is non-uniform, the usual definition of rg needs to be further clarified 
as to what value of p should be used. In our case, it is convenient to set p = p out — p(r = r out ), i.e. we use the gas 
volume density at the outer cloud core boundary. Then, the resulting free-fall time of the outermost gas layer becomes 

Tff.out = V37r/(32Gp out ). 

We now need to relate rg iOUt with the initial cloud core mass M c \, Neglecting a small plateau in the gas surface 
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density profile at r < r Q , we can write (see Section [3]) 



£ = 



The mass of a cloud core with this gas surface density distribution is 

2A 1 /2 c 2 rout 2A4 



(Al) 



(A2) 



G 7rG£ ut 

where S out is the gas surface density at the cloud core's outer boundary. Assuming a local vertical hydrostatic equi- 
librium with the scale height Z = c 2 /(irGE), approximating the gas volume density as p = T,/(2Z), and substituting 
the resulting value of p ou t = T, 2 ut irG / (2c 2 ) in the expression for Tff iOUt , we finally obtain 

T ff ,out[Myr] = ^(^! "l^l, (A3) 



JLY' ( Mcl \ 

Jo) \M Q J ' 

where To is the initial gas temperature, p is the mean molecular weight, and A is the amplitude of initial density 
enhancement. 



